Classical Strongly Coupled QGP I : 
The Model and Molecular Dynamics Simulations 
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We propose a model for the description of strongly interacting quarks and gluon quasiparticles at 
T = (1 — 3)T C , as a classical and nonrelativistic colored Coulomb gas. The sign and strength of the 
inter-particle interactions are fixed by the scalar product of their classical color vectors subject to 
Wong's equations. The model displays a number of phases as the Coulomb coupling is increased 
ranging from a gas, to a liquid, to a crystal with antiferromagnetic-like color ordering. We analyze 
the model using Molecular Dynamics (MD) simulations and discuss the density-density correlator in 
real time. We extract pertinent decorrelation times, diffusion and viscosity constants for all phases. 
The classical results when extrapolated to the sQGP suggest that the phase is liquid-like, with a 
diffusion constant D » 0.1/T and a bulk viscosity to entropy density ratio rj/s ~ 1/3. 



I. INTRODUCTION 

The Quark-gluon plasma (QGP) is a high tempera- 
ture phase of QCD. Here, the word plasma is used in 
the same sense as in electrodynamical plasma, that is a 
phase with free color charges that are screened [1] rather 
than confined by the medium. Lattice simulations have 
shown, that the QGP exists above the critical tempera- 
ture T > T c « 170 MeV of the QCD phase transition, in 
the deconfined and chirally symmetric phase. 

Asymptotic freedom of Non-Abelian gauge theories, a 
property of QCD, insures that for high enough temper- 
ature T A.qcd this phase becomes weakly coupled 
(wQGP), with most of the particle interactions charac- 
terized by a small coupling a s {p w T) -C 1. In this 
regime, the wQGP is a near-ideal gas of its fundamen- 
tal constituents, quarks and gluons. Perturbative meth- 
ods, such as hard thermal loops [2] and other resumma- 
tions techniques, show how quarks and gluons are dressed 
and become quasiparticles with T-dependent dispersion 
curves and widths. 

The QGP is experimentally studied using heavy ion 
collisions, such as dedicated experiments carried at the 
CERN SPS during the previous decade, and now at the 
BNL RHIC collider. RHIC can reach temperatures of 
about 2T C . The success of the hydrodynamical descrip- 
tion for the observed collective flows at RHIC [3] has 
shown that all dissipative lengths are very short. The 
produced matter at RHIC cannot be a weakly coupled 
gas but a rather good liquid [4]. Recently, two of us [5] 
have suggested that the interaction of quasiparticles in 
the relevant temperature range at RHIC is strong enough 
to generate multiple marginal colored Coulomb bound 
states. Some of those states (charmonium) are observed 
in current lattice [6] simulations till about 2.5T C . The ex- 
istence of these states, especially the colored ones, is still 
debated and warrant further numerical and independent 
checks on the lattice. 

The effective potential energy of two static colored 
charges separated by a distance R, U(R) can also be 
deduced from lattice simulations. Close to the crit- 



ical temperature, the separation energy AU (the po- 
tential at infinity minus its value at some typical dis- 
tance 0.3 fm) is AU w 4 — 1 GeV in the temperature 
range T = (1 — 1.2) T c = 0.17 - 0.21 GeV. The ensu- 
ing Boltzmann penalty e~ u with u — AU/T w 20 — 5 
is mighty, indicating the dominance of potential over ki- 
netic energy. This regime is now called a strongly coupled 
QGP (sQGP). Its structure and consequently its trans- 
port properties are radically different from a wQGP to 
which it is expected to fold at very high temperature. 

In the traditional context of an electromagnetic 
plasma, the term "strongly coupled" plasma has a similar 
meaning. Ionic or dusty plasmas have charged ions with 
large masses, and thus are essentially classical. The stan- 
dard dimensionless parameter characterizing the strength 
of the interparticle interaction in a classical plasma is T 
the ratio of potential to kinetic energy 
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where Ze,aws>T are respectively the ion charge, the 
Wigner-Seitz radius ay/s = (3/47rn) 1 / 3 and the tempera- 
ture * . r is convenient to use because it only involves the 
input parameters, such as the temperature and density. 
However, one should keep in mind that the real interac- 
tion parameter is the output dimensional parameter u, 
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T 



(2) 



where U =< V > is the average interaction of a particle 
with all its neighbors. This is the ratio that enters the 
effective Boltzmann exponent, and defines all correlation 
functions. 

Since u is proportional to T, one usually defines the 
weakly coupled regime for r < 1 and the strongly cou- 
pled regime in the opposite limit f> 1. Extensive stud- 



*We do not use fcs and thus measure temperature in energy 
units. 
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ies of the one-component plasma (OCP) in electrodynam- 
ics, using both MD and analytical methods over the past 
decades, have revealed the following regimes: i. a gas 
regime for T < 1; ii. a liquid regime for T w 10; iii. a 
glass regime for T w 100; iv. a solid regime for T > 300. 
For a review see e.g. [8]. 

Another physical system closer to what we will discuss 
in this work is a classical two-component plasma (TCP) 
with both positive and negative charges. Examples are 
molten or frozen salts with ions of comparable masses, 
as in the sQGP. Hydrogen plasmas are better studied, 
but the underlying charges carry very different masses, 
so the ensuing results bear no insight to the sQGP were 
the masses are likely comparable. 

Quantum effects will be schematically incorporated 
in the form of a localization energy, providing a repul- 
sive core for the 2-particle interaction irrespective of the 
charge. In the future we plan to study more quantum 
corrections to classical MD, in particular the role of the 
exclusion principle, and the relation between classically 
correlated charges and quantum bound states. 

We are certainly aware of many cases of quantum plas- 
mas. An quantum electron plasma in metals is quan- 
tum and thoroughly studied, but repulsive nature of the 
electron-electron interaction excludes the formation of 
bound states and is very different from QGP we are in- 
terested in. Only weak induced electron-electron interac- 
tion leads Cooper pairing in superconductors. Excitons 
in semiconductors, in a plasma of particles and holes, are 
even closer to our problem. Recent progress in trapped 
strongly coupled fcrmionic atoms have further elucidated 
how the transition to Bose-condensed pairs takes place. 
However, in all these cases one needs novel tools to study 
quantum effects dynamically: the existing ones, such as 
e.g. restricted paths Monte Carlo, use Euclidean time 
correlators which (like those from lattice QCD) are good 
for thermodynamical observables but next-to-impossiblc 
to use for transport properties. The MD approach, what- 
ever crude it can be, is an invaluable classical tool which 
directly provides real-time correlators. Furthermore, per- 
tinent quantum corrections can be added later, as it is 
done in some atomic systems. 

The model we propose in this work combines features 
of strongly coupled Abelian TCP, with non-Abelian fea- 
tures proper to QCD in the form of classical color vectors 
for the underlying charges. The main interparticle inter- 
action is proportional to the dot-product of the color vec- 
tors, which of course can be attractive, repulsive or even 
null if the color vectors of the particles are orthogonal to 
each other. Dynamics of color vectors as well as particle 
coordinates are described by classical equations of mo- 
tion (EoM), and thus can be studied by MD simulations, 
which are less involved than the corresponding coupled 
hierarchy of quantum Green's functions. Below, we will 
argue for which values of the parameters this model can 
be useful for understanding the sQGP. Before we will do 
so, the model will be defined and studied in its own right. 



II. CLASSICAL QUARK-GLUON PLASMA, 
CQGP 

A. The model 

The model is based on the following main assumptions: 

i. The particles are heavy enough to move nonrelativis- 
tically, with masses M » T; 

ii. The inter-particle interaction is dominated by colored 
electric (Coulomb) interactions, with all magnetic effects 
(like e.g. spin forces) ignored; 

iii. The color representations are large, so that color op- 
erators t a can be represented by their average, classical 
color vectors. 

The parameters of the model include the particle mass 
M and their density n. The main interaction potential is 
proportional to the dot-product of the unit color vectors 
times the standard Coulomb interaction strength. The 
notations we use parallel the conventions in EM plasmas 
and in QCD. The relation between the two is defined as 

(Z a ef=C a ^- (3) 

where Z 2 in EM is C a in QCD, the Casimir operator 
eigenvalue for q,q,g, while e 2 = g 2 /47r t. When needed, 
we will "input" the Debye screening mass" Md, if the 
classical screening of plasma charges derived dynamically 
from MD would not be an accurate representation of the 
effective forces between quasiparticles in the sQGP. 

In QCD, all the effective parameters used here are func- 
tions of the temperature T and the baryon chemical po- 
tential ji. or of a point at the QCD phase diagram. In 
a heavy ion collision, those are defined by the collision 
energy and centrality of the collision. Furthermore, for 
each volume element the cooling of the matter during its 
expansion is well approximated by the adiabatic (fixed 
entropy/baryon densities) path on the phase diagram, 
related T and /U. 

As a first approximation, baryonic charge play little 
role at RHIC and thus one can set \x = and think of 
only one parameter, the temperature T, defining matter 
properties. Therefore, when we attempt to map the MD 
results for the cQGP into sQGP, all the parameters of 
the model mentioned above should be converted to tem- 
perature, as we will detail. In way, the classical model 
we are about to discuss in an unrestricted n-dimensional 
parameter space, will be useful for the sQGP through a 
1-dimcnsional slicing of the n-dimensional space. Prob- 
ing the n-dimcnsional parameter space for the model is 
useful theoretically. Much like EM plasmas, the dimcn- 
sionless parameter (1) will be key in characterizing the 



TNote that the standard notations used in EM and QCD 
differ by 4tv in the Lagrangian and all subsequent formulae. 
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transition from a weakly coupled regime with r -C 1, to 
a strongly coupled regime r>l. 

At very low temperatures T or very large T, any classi- 
cal system will freeze. Therefore, the low temperature be- 
havior is dominated by the lowest ordered state. Here, we 
recall that in general, the Abelian 2-component plasma 
freezes to an ionic crystal much like ordinary salt NaCl. 
The non- Abelian plasma under consideration freezes also 
to cubic crystal with ferromagnetic (alternating) order of 
the classical color vectors. A direction in classical color 
space is selected randomly, through spontaneous symme- 
try breaking of the global color group, with alternating 
directions of the color vectors along the crystal axes. 

B. The equations of motion 

The specific model Hamiltonian is 

H = E IT*- + Vc + V «™ ( 4 ) 

a i 

with the standard kinetic energy, the colored Coulomb 
interaction Vc and the repulsive core V COTe respectively. 
The colored Coulomb part is 

a i=tf3 i ^ J 

The sums are over the species a = q,q,g and their re- 
spective numbers N a . We will specify V corc below. 

Thus our phase space coordinates are the position (x a ), 
momentum (p a ) and color (Q a ). The latter rotates under 
a gauge transformation 

Qa ^ D ab^ Q b (g) 

showing that (5) is gauge invariant. Only the Coulomb- 
like interaction was retained in (5) since the magnetically 
induced interactions are subleading non-relativistically. 
This is the case for all spin and local many-body forces, 
for instance the local 3-body gluon-induced interaction 

E Q a at Q%Q c lk ( f abc F + d abc r>) (7) 

is subleading F w D w V/m. Non-local many-body 
interactions as induced by the 2-body Coulomb and core 
interactions are important and will be resummed to all 
orders using MD. 

The EoM for the phase space coordinates follow from 
the usual Poisson brackets. For the standard coordinates 
they are 

{x^ j } = S mn 5 afi 5 ij (8) 
For the color coordinates they are 

rc„^Hf k Q« I (9) 



the classical analogue of the SU(N C ) color commutators, 
with f abc the known structure constants of the color 
group. The classical color vectors are all adjoint vec- 
tors with a = l...(A 2 — 1). For simplicity only the non- 
Abelian group SU(2) will be considered, with 3d color 
vectors. The difference between quarks and gluons fol- 
low from their respective Casimir assignments. 

Although the brackets (9) do not look like the usual 
canonical relations between coordinates and momenta, 
they actually are. Indeed, a more standard phase space 
description requires the use of the Darboux parameteri- 
zation [9] 

Q=($ A ,IL A ) (10) 

with A = 1,..,N C (N C — l)/2 now satisfying standard 
canonical commutation relations 

{<f> A ,n B } = 5 AB (11) 

the $'s are angular coordinates and the IT's are angular 
momenta. The ll's are identified with the fixed Casimirs 
of SU(N C ) which are conserved by the equations of mo- 
tion. For example, for N c = 2, the 3 components of the 
color vector conserve one Casimir operator C2 (the color 
vector's length). Without loss of generality, the Darboux 
coordinates maybe chosen as 

(0 1 ±zQ 2 ,Q 3 ) = (v / J 2 -n 2 e ±4 *,n) (12) 

with J 2 /3 = C 2 . For N c = 3 there the color vec- 
tor is 8-dimcnsional, with 2 fixed Casimirs (Q a Q a ) and 
(d abc Q a Q b Q c ) 1 and 6 variables. In this case, the Darboux 
set is more involved, but can be parametrized with 3 an- 
gles and 3 conjugate momenta. We will not be working 
with Darboux-based EoM below. 

The pertinent equations of motion for our classical 
Coulomb gas resulting from (4) and the Poisson brackets 
are 

Pa i 

m a 

Q a a * = {H,Q al } = g.f abc Q b at A c at0 (13) 

with 

fpa v7 A a 

i — v «' i 

The last of the three equation of motion in (13) is also 
known as Wong's equation [10]. Due to the antisym- 
metric nature of the structure constant f abc , it has the 
form of rotation of the color vector conserving its length 
(classical precession). For SU(2) this is the only Casimir, 
while for higher color representations more Casimirs are 
involved. We expect the resulting EoM to generate a bal- 
let of classical precession conserving all these Casimirs. 
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III. NUMERICAL STUDIES USING 
MOLECULAR DYNAMICS 

In this paper we will outline the basic ingredients in- 
volved in the MD analysis of (13). We first discuss what 
the pertinent scales are and how they are set. Second, we 
briefly comment on the numerical analysis used to deal 
with long-range Coulomb forces. 



the quantum localization is the introduction of a dimen- 
sionful parameter, i.e. hbar in natural classical units. 

i. Length unit : At any n, the minimum of the 
potential is at r = A. So A sets the basic length scale, and 
from now on all distances are measured in units of A. For 
example, a close cubic-packing density with n cp = 1/A 3 
will be written as n cp = 1. ii. Time unit : The unit of 
time ro will be set by the natural frequency of a screened 
Coulomb system, its plasma frequency 



A. Potential and units 

The MD simulation is basically a numerical solution of 
EoM starting from some initial conditions. The method 
was developed in the early days of computers in the 
1950's, and applied extensively for the studies of sim- 
ple liquids and plasmas. It is the most straightforward 
way to access transport properties such as the diffusion 
coefficient D, the bulk viscosity rj, the heat conductivity 
A, as well as pertinent decorrelation functions and times. 

In standard 2-component plasma physics MD simula- 
tions, the interparticle potential is chosen as 



TO = LU V 



V(r) = V COIC + V C = 



-(- 

n \ r 



(15) 



with r = |a?j— x*j\ and Q i7 Qj — ±1 are respective charges. 
The repulsive core potential V core ensures overall mat- 
ter stability with n some large parametric number. For 
atomic systems n — 9. For molten salts the core follows 
from the intrinsic repulsion of the atomic electrons. 

Since this is the first paper in a series on this model, 
and as we arc mostly interested in understanding the 
role played by the non-Abclian color variables, we have 
carried the numerical MD analysis for both the Abelian 
and non-Abelian charges. Both analysis were carried in 
parallel to understand the differences between discrete ± 
charged plasmas and our classically colored SU(2) plasma 
with Q a charges. Therefore, the same Vcorc was kept in 
the non-Abelian simulation. 

The microscopic motivation in QCD for V'coro stems 
from the short range part of the inter-particle interac- 
tions. Needless to say that this is a hard question, and 
not much is known about it from first principles whether 
in the intermediate coupling region (a s w 1/2) or the 
very strong coupling region (a s k 1 or T > 1). 

One reason for an effective repulsion at short distance 
is the so called quantum localization potential 



mr* 



(16) 



Such potentials are usually used in classical studies/ sim- 
ulations of light atoms, e.g. solid He. Note that this term 
has a lower power n = 2 rather than n = 9 used in V COIC - 
This not a much importance as any n > 1 stabilizes the 
Coulomb attraction at short distances. The net effect of 



Airne 2 



1/2 



(17) 



iii. Mass unit : This is naturally given by the parti- 
cle mass m. Thus all important time correlators have a 
characteristic length of order 1 in To. In contrast, all MD 
runs are typically few hundreds or thousands tq to insure 
statistical equilibrium. 

This fixes the three basic units used in our analysis. 
For instance, the kinetic energy is measured in obvious 
units of m\ 2 /tq. The strength of the Coulomb potential, 
if measured in such units, is 



A Tq 



1 1 

47T 7lA 3 



(18) 



so that it is defined directly by the particle density. The 
temperature, in the same units, can be defined via av- 
eraged (dimcnsionless) velocities as (3/2)T =< v 2 /2 >. 
Finally the main parameter u in (2) can be obtained from 
simulation, as the ratio of the measured total mean po- 
tential to the total kinetic energy. 



B. Technical details 

There are readily available numerical packages for solv- 
ing the set coupled set of differential equations of the 
type (13). In our case, their accuracy is tested through 
energy conservation and color vectors remaining on the 
unit sphere. We have used the package from the CERN 
Library in double precision and tested its accuracy for 
both energy and color conservation. 

The many particle problem in a box requires the in- 
troduction of boundary conditions. We chose periodic 
boundary conditions. When a particle crosses one side, it 
reappears on the opposite side. We can visualize space as 
filled with many mirror cubes, piled next to each other. 
When one particle crosses out of the central cube, its 
mirror image enters from the opposite cube. The inter- 
particle potential in this many-cube world is periodic un- 
der cubic translation by L, i.e. 



^ V(Si — xj + n L) 



(19) 



n x ,n v ,n z 



So as a particle crosses the boundary, there is no change 
in the force or the potential since (19) is truly periodic. 
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In practice, the sum in (19) is only carried over few ad- 
jacent cubes or mirror images, causing an overall loss of 
periodicity. As a result, the forces at the faces of the 
cube jump slightly thereby causing random kicks on the 
crossing particles. This gradually heats the system and 
should be controlled. The more long range the potential 
is, the more severe the heating. 

Traditionally the Bethe polynomials for Ewald sums 
are used to overcome this heating phenomenon caused 
by the truncation. We have found this procedure not 
sufficiently accurate for the amount of complexity it in- 
troduces. Indeed, it leads to a factor of 50 suppression of 
the electric field normal to the cubic interface. Current 
computers easily allow for keeping up to 3-4 mirror im- 
ages on each side of a central cube, thereby reducing the 
net electric field by orders of magnitude. 

Another practical way to deal with this problem is in- 
troduction of a very small friction into the EoM, which 
can compensate the heating effect caused by the random 
forces at the cubic interface. Our typical runs consists of 
4 4 = 64 or 6 4 = 1296 particles all of the same species. 

C. Structure factor 

The static and dynamic properties of a liquid are deter- 
mined by various correlations functions. One of the most 
useful of these functions is a density-density correlation 
function or structure factor defined as, 

G(r,t) = ±(p(x,t)p(0,0)) , (20) 

with r — \x\ and p(x,t) being the particle density at the 
position x at time i, 

JV 

p(f 7 t) = Y / S(x-x t (t)), (21) 

i=l 

Here, n = N/V is the particle density and the averaging 
in (22) is carried over an equilibrium ensemble *. 

Assuming translational invariance and isotropy, the 
classical density-density correlation function can be writ- 
ten as, 

/ N N V 

G ^ = ]v ( ^ 15 (f + f l(0) ~ Sj {t)) / ' (22) 

with TV is the number of particles, Xi(t) is the position 
of the ith-particle at time t. G(r, t) characterizes the 
likelihood to find 2 particles a distance r away from each 
other at time t. 
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FIG. 1. Gd correlation function for T = 0.83,31.3, 131, re- 
spectively. Red circles correspond to t* = 0, and blue squares 
correspond to t* = 6. 



*In a given MD simulation the total energy and number of 
particles of a system is kept constant. This corresponds to 
using a microcanonical ensemble. 



The density-density correlation function (22) can be 
written as a sum of two parts referred to as the van Hove 
or self- correlation G s , and the distinct- correlation Gd, 

G(x,t) = G s (x,t) + G d (x,t) , (23) 
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where 



G s (x,t) = jj(fl S(x + Xi(0)-Xi(t))J , (24) 
and, 



We have measured a large selection of such function, 
as a typical examples let us present here Gd at three 
values of the coupling, moderate , strong and very strong, 
r = 0.83,31.3,131, see Fig.l. One sees from those that 
in the first case a relatively weak correlation between the 
particles, at distance 1 corresponding to the potential 
minimum, which relaxes rather quickly with time. The 
correlation is more robust in the second "liquid" case, 
and is very stable and is accompanied by extra peaks in 
a "solid" last case. 



D. Transport coefficients 



An important aspect of the strongly coupled plasmas is 
their dramatic change in transport properties in compar- 
ison to weakly coupled plasmas. The current MD simula- 
tions can be used to study the bulk transport properties 
at strong coupling. In particular, the shear and bulk vis- 
cosities, diffusion constant, thermal conductivity, color 
conductivity etc. A simple way to obtain these transport 
coefficients is via Green-Kubo relations. These relations 
give the transport coefficients in terms of the integrals of 
the equilibrium time-dependent correlation functions. 

In the case of the self-diffusion the corresponding cor- 
relation function is the velocity autocorrelation function, 

where Vi(r) is the velocity of a particle i at time r. The 
velocity autocorrelation functions are shown in Fig. 2 for 
r = 0.83,31.3,131 respectively. This refers to the non- 
ideal gas, liquid and crystal of the (one species) cQGP. 

The diffusion constant is the integral of the velocity 
autocorrelation function, 

/>oo 

D= / D(r)dT. (27) 
Jo 

Fig. (3) shows a log-log plot of the diffusion constant 
as a function of T . The dependence of D on T is linear 
and can be approximately described by a simple power 



(28) 
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FIG. 3. The diffusion constant of a one species cQGP as a 
function of the dimensionless coupling V. Blue points are the 
MD simulations; the red curve is the expression (28). 



The viscosity coefficients are given in terms of time 
autocorrelation function of the stress-energy tensor, 



1 



3TV 



(29) 



\x<y 



where V denotes a sum over the three pairs of distinct 
tensor components (xy,yz and zx). The off-diagonal 
parts stress-energy tensor are given by 



N I 



(30) 



and cyclically, with fi 



Xj and Fi-i is the force 



on particle i due to particle j. The stress-tensor correla- 
tion functions are shown in Fig. 4 for T = 0.83, 31.3, 131 
respectively. 

The Green-Kubo relation for the the coefficient of the 
shear viscosity is, 



V 



t)(t) dr . 



(31) 



The coefficient of shear viscosity as a function of T is 
shown in fig. 5. For small V the viscosity is large since 
the mean-free path is large in a gas-like phase. The vis- 
cosity is minimum in the liquid phase, and rises slowly in 
the crystal phase at large V. The coefficient ij is approx- 
imately fitted by 





FIG. 4. Stress-tensor autocorrelation function t)(t) for 
T = 0.83, 31.3, 131. 



ti re o.ooi r 



0.242 0.072 



p0.3 



r 2 



(32) 



Notice that the fit misses 77 « 0.1 at T = 10 which is the 
lowest viscosity measured in the current cQGP. 
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FIG. 5. Shear viscosity of one species cQGP as a function 
of the dimensionless coupling F. Blue points are the MD 
simulations; the red curve is the fit in eq. (32). 



while that of the decorrelation time is 
t v w 0.239 + 0.091 VT. 



(35) 
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The stress-tensor auto-correlation function can be used 
to determine the viscous decorrelation time t, ; , defined as 



V 



'7(0) 



(33) 



where 77(0) is the value of the stress-energy auto- 
correlation function at time t = 0. The values of 77(0) 
and Trj as a function of T are shown in Fig. 6 and Fig. 7, 
respectively. 



T|0' 




FIG. 6. Values of stress-energy auto-correlation at zero 
time of a one species cQGP as a function of the dimensionless 
coupling F. The blue points are the MD simulations; the red 
curve is the fit in eq. (34). 



The functional dependence on T of 7/(0) is 

0.77 0.44 



r?(0) « 0.0005 r 



11.57 



po.33 



(34) 



FIG. 7. The viscous decorrelation time 
for the stress-energy tensor time correlator as a function of 
the dimensionless coupling F. The blue points are the MD 
simulations; the red curve is the fit in eq. (35). 



E. Potential energy 

Another useful observable that can be determined from 
MD simulations is the total potential energy U normal- 
ized to T. This ratio as a function of T is shown in Fig. 8. 
The parametric dependence of this ratio on T is given by, 

U u 2.2 (36) 



-4.9 - 2r + 3.2r 1/4 



nt " " " ' r x /4 ■ 

At large T the crystal phases sets in, and the potential 
energy asymptotes U '/NT » — 2r which is a measure of 
the Madelung constant. A full theoretical analysis of (36) 
will be given in the sequel II. 
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FIG. 8. Ratio of the potential energy to NT of a one 
species cQGP as a function of the dimensionless coupling Y. 
The blue points are the MD simulations; the red curve is the 
fit in eq. (36). 



IV. COMPARISON WITH SQGP 

A qualitative mapping of the one species cQGP de- 
scribed above into sQGP can be made with the by ad- 
justing the three basic scales described above, namely 
that of the length, the time and the mass. As we have 
repeatedly indicated above, all parameters of the model 
are functions of the temperature T in the sQGP. In this 
section we provide heuristic arguments for what the pa- 
rameters are using a simple description of the sQGP at 
T = 1.5 — 3 T c . (It is used as a reference point because 
more lattice data are available for it than for T close to 
T c .) 

As discussed above, the unit of length is A the mini- 
mum of the potential. Since the repulsive part mocks up 
the quantum repulsion at short distances, the effective 
inter-particle potential is § 



2mr 2 



(37) 



with a minimum at r = K /(mCa s ) = A. C is the 
pertinent Casimir for quarks and gluons. 

The quasiparticle mass, defined as its energy at zero 
momentum, is at T = (1.5— 3)T C about constant to 5T C 
both for quark and gluon quasiparticles [12]. We note 
that there are no direct lattice measurements of quasi- 
particle dispersion curves below 1.5 T c . One can infer 
quark quasiparticle masses from fits to thermodynamical 
observables and baryonic susceptibilities: it tells us that 
masses grow toward T c , for a recent discussion see [13]. 

For further discussion it is convenient to introduce a 
dimensionless mass parameter m = m/T, which thus 
changes from large values close to Tc to about 3 at 1.5T C . 
Clearly, its value determines how accurate is the nonrcla- 
tivistic approximation. For example, at to = 5 the mean 
square momentum \f< p 2 >/T = 4.8 if calculated rela- 
tivistically, while the non-relativistic approximation gives 
the value smaller by a factor 0.8. The corresponding ra- 
tio of the total density in non-relativistic approximation 
is .7 of its relativistic value. So, it is not an unreasonable 
approximation, but not very accurate. At T = 1.5T C the 



§ In so far, our arguments were classical and nonrelativistic, 
so both h and 1/c were set to zero. In this section, and 
this section only, we use the standard high energy units with 
h = l,c = 1. Hopefully, it will not create confusion for the 
reader. We will keep Ti in the next formula, to emphasize the 
quantum origin of the repulsive core. 



value of this mass is m w 3 only, and the accuracy of it is 
even worse. Therefore, we should calculate any momen- 
tum integrals relativisitcally, and only the results to be 
mapped into cQGP. 

The color Casimir for gluons C g = 3 and quarks C q = 
4/3 should be averaged over their respective weights. The 
number of effective degrees of freedom in thermodynam- 
ical quantities is such that roughly all three species - 
g,q,q - are equally represented. Furthermore, the values 
of the coupling constant a s inferred from measured static 
potentials at relevant distances is a s rts 0.5. Withing the 
uncertainties, we will thus simply use < a s C >= 1 be- 
low. Thus 



A = r 



3T 



(38) 



The density of quasiparticles n can be estimated as 
follows. Although quasiparticles are relatively heavy, the 
presence of bound states compensates and increases the 
density of total quasiparticles to about 0.8 of the entropy 
density per particle in the massless gas [14]. Ignoring this 
0.8 which refers to the bound states, we then have 



(0.244 T 3 ) (8 + 6 TV/) w 6.3 T 3 



(39) 



where the first bracket is the usual density of black body 
radiation photons, and the second bracket is the number 
of effective degrees of freedom due to 8 colored gluons and 
2 quarks and antiquarks with 3 colors and Nf flavors This 
corresponds to the following Wigner-Seitz radius 



a-ws 



3 



4 irn 



1/3 



1 

3T 



(40) 



The time units in MD are given by the plasma frequency, 
which is 



T = UJ V 



Ann 



a s C> 



TO 



1 



5.1T 



(41) 



where, for definiteness, we used m = 3T. 

In summary, the dimensionless cQGP results of the 
MD carried above can be qualitatively translated to the 
dimcnsionfull sQGP parameters using respectively the 
length A, time tq and mass to units 

As " Jt 
i 



Tn ~ 

5.1 T 

m w 3T 



(42) 



with r w 3. 

Indeed, the viscosity unit rj is then given by the com- 
bination of units with the appropriate dimension 

r?o = mj (r A) - (3T)(5.1T)(3T) « 46 T 3 (43) 

leading to the sQGP viscosity through 
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7] -> V r] w (0.17)(46T 3 ) w 7.8 T 3 



(44) 



In the sQGP, it is customary to give the viscosity per 
entropy density, which for ideal massless QGP is 



x 2 x 2 x 3 x N f T 3 w 20T 3 (45) 



for 3 flavors. Correcting a bit for more realistic lattice 
entropy, we will use s w 23 T 3 around T = l.bTc. Thus, 
the dimcnsionless viscosity r\ w 0.17 as measured in the 
one species cQGP above at T rj 3, can be translated 
to the dimensionfull viscosity in the sQGP by rescaling 
through t] . The viscosity per entropy ratio in the sQGP 



is 



1 
s 



Wo 
s 



0.34 



(46) 



N=4 SYM results have suggested a universal lower bound 
for rj/s > 1/Att. Our analysis suggests that in the sQGP 
this ratio is about 4 times the lower bound. 
The diffusion constant translates to 



D w 0.161 -» (0.161)(A 2 /t ) « 0.1/T . 



(47) 



Similarly, the stress-tensor decorrelation time translates 
to 



a liquid-like phase at intermediate coupling, and fi- 
nally a crystal-like phase at strong coupling with anti- 
ferromagnetic-like color ordering. The transition from 
liquid to crystal is best seen in the density-density cor- 
relation function or its static counterpart the structure 
function. At large Coulomb coupling the excess energy 
per particle asymptotes the Madelung energy for a crys- 
tal. 

We have used the numerical analysis to extract a num- 
ber of transport coefficients including the quasiparticle 
diffusion constant, the bulk viscosity and relaxation time. 
While quantum mechanics is important in the sQGP in 
the formation of the underlying quasiparticles (dispersion 
laws) and the running of the coupling, we have suggested 
that the ensuing quasiparticle interaction and dynamics 
is essentially Coulombic and at strong coupling classical. 
For this, the cQGP results are generic at intermediate 
and strong coupling. We have given qualitative argu- 
ments for how to relate the cQGP transport results to 
those of interest in the sQGP by identifying the pertinent 
length scales. Indeed, we have found that the sQGP cor- 
responds to T w 3 which is liquid-like and viscous with 
rj/s w 1/3T. 

In a series of sequels to follow, we will provide more 
insights to the cQGP results presented above for the bulk, 
2-particle correlations and transport properties. 



t v w 0.395 -> (0.395)(t ) « 0.08/T 



(48) 



One should keep in mind, that these results are estimated 
for T s=s 1.5T C . As temperature gets closer to T c , the 
masses of quasiparticles grow further and the coupling 
gets larger, so one finds even stronger coupled plasma. 
Unfortunately, we do not yet have sufficiently precise lat- 
tice data to make the mapping of the cQGP to the sQGP 
in this region more definite. 
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V. CONCLUSIONS AND PROSPECTS 

Quarks and gluons in the temperature range of 1-1.5 T c 
behave as quasiparticles with masses larger than m > 3 T 
and a Casimir-weighted Coulomb strength of the order 
of < C a s 1 which is strong. Lattice measured po- 
tentials in this temperature range give AU /T between 
20 at T c and 5 at 1.2T C indicating the dominance of the 
potential energy over the thermal energy. 

In this paper, we have suggested that the long- 
wavelength properties of the QGP in the 1-1.5 T c range 
may be modeled by a classical and non-relativistic gas of 
massive quasi-particles interacting via strong but classi- 
cal color charges. The model we coined cQGP has merits 
on its own and was analyzed numerically using molecu- 
lar dynamics together with Wong's equations for classical 
color evolution. 

Our results show the existence of several phases rang- 
ing from a gas-like phase at weak coupling, through 
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